/*********************************************/
/*                                           */
/*  This program does post-estimation        */
/*  analysis of the worker and firm          */
/*  fixed effects estimated for our          */
/*  AKM/CHK analysis.                        */
/*                                           */
/*  Before running this we need to make      */
/*  the appropriate text file printouts      */
/*  of the AKM/CHK estimation. These         */
/*  printouts include all the estimated      */
/*  fixed effects                            */
/*                                           */
/*  Also need a file that has firm size      */
/*  data to merge in here.                   */
/*                                           */
/*                                           */
/*********************************************/
%let basedir = "<basedir>"; /* matlab estimation output dataset */
%let outputdir = "<outputdir>";  /* where the output files go        */
%let firmdir = "<firmdir>";  /* where the firm files are stored   */
%let sample = S100P;           /* Sample is either 10% (S10P) or 100% (S100P)  */
%let spec = 1;               /* =1 if earnval=520hr */
%let labsuprest = 0;         /* 0 = no labor supply restriction, 1 = sandwich */
%let earnval = 520hr;       /* Minimum number of hours worked , for sandwich "520hr_Sandwich" */
%let genderval = Male;           /* Population is either: Male, Female, AllSex */
%let minwageval = DEFLATE;           /* Minmum wage restriction is either: CTEMP (contemporaneous) or DEFLATE (2013 PCE deflated) */

/**********************************************************************/
/*                                                                    */
/*  COMPUTE POST-ESTIMATION STATISTICS                                */
/*                                                                    */
/*                                                                    */
/*    There are no firm size, industry, or job switcher per firm      */
/*	  restrictions. Instead we compute all the relevant moments       */
/* 	  within these groups so that any aggregate relationship can      */
/*	  be computed with the conditional moments and # of cell obs.     */
/*                                                                    */
/**********************************************************************/
libname firmdata "&firmdir.";
%macro feanalysis(earn,gender,minwage);
%do K=1 %to 5;   
%if &K=1 %then %do;%let FY=1980; %let LY=1986;%end;
%if &K=2 %then %do;%let FY=1987; %let LY=1993;%end;
%if &K=3 %then %do;%let FY=1994; %let LY=2000;%end;
%if &K=4 %then %do;%let FY=2001; %let LY=2007;%end;
%if &K=5 %then %do;%let FY=2007; %let LY=2013;%end;   /* program runs for all 5 estimation periods */
/**********************************************************************/
/*                                                                    */
/*  0.) a.) Import AKM output data.                                   */
/*      b.) Merge in firm data.                                       */
/*      c.) Create FE decile, industry, and firm size groups           */
/*                                                                    */
/**********************************************************************/
proc import datafile="&basedir.\&sample._AKMeffs_&spec._&labsuprest._&minwage._&earn._&gender._Int7_&K..txt" /* Import fixed effect output. Check directory above. */
	out=fe_out
	DBMS=TAB
	replace;
	getnames=no;
run;
data fe_out;                   /* label variables */
	set fe_out;
	rename var1=id var2=firmid var3=year var4=pe var5=fe var6=xb var7=r var8=y var9=yob var10=ind var11=geo var12=n_jobs var13=i_sw;
	age=var3-var9;
run;
proc sort data = fe_out;by year firmid; run;
%macro import;
%do J=&FY %to  &LY;  
data firmsizes&J.;                       /* import firm size data */
	set firmdata.firmsizes&J.;
	year=&J; 
	keep year ein_original countppl indcode geo;
	rename ein_original=firmid;
	rename countppl=firm_size;
	rename indcode=ind;
run;
%end;
%mend;
%import;
data firmsizes;
	set %do J=&FY %to &LY; firmsizes&J %end;;
run;
proc datasets lib=work nolist;
	delete %do J=&FY %to &LY; firmsizes&J %end;;
quit;
proc sort data=firmsizes; by year firmid;
data fe_out;                       /* merge in firm data */
	merge fe_out (in=p) firmsizes;
	by year firmid;
	if p=1;
run;
/****************************************/
/*  Compute the AKM match effect        */
/****************************************/
/*
proc sort data=fe_out; by id firmid; run;
proc means data=fe_out noprint nway;
	class id firmid;
	var r;
	output out=match mean(r)=m;
run;
*/
/* merge back to main data */
/*
data fe_out;			
	retain id firmid year yob age y pe fe xb r m e;
	merge fe_out match;
	by id firmid;
	drop _freq_ _type_;
	e = r-m;
run;
proc datasets lib=work nolist;
	delete match;
quit;
*/
/***************************************/
/*  Compute the number of periods      */
/*  each worker is in the interval     */
/***************************************/
/*
proc sort data=fe_out; by id year; run;
data n_periods;
	set fe_out;
	by id;
	x=1;
	if first.id then n_wper=0; n_wper+x;
	if last.id;
	keep id n_wper;
run;
*/
/* merge back to main data */
/*
data fe_out;   
	merge fe_out n_periods;
	by id;
run;
proc datasets lib=work nolist;
	delete n_periods;
quit;
*/
/***************************************/
/*  Condense industry codes to         */
/*  12 categories                      */
/***************************************/
data fe_out;
	set fe_out;
	if ind^=0 then ind=floor(ind/100);
run;
data fe_out;
	set fe_out; /* SIC codes */
	if 1<=ind<=9 then ind=1; /* Agriculture, Forestry & Fishing */
	if 10<=ind<=14 then ind=2;  /* Mining */
	if 15<=ind<=17 then ind=3;  /* Construction */
	if 20<=ind<=39 then ind=4;  /* Manufacturing */
	if 40<=ind<=49 then ind=5;  /* Transportation, Communications, Electric, Gas & Sanitary Services */
	if 50<=ind<=51 then ind=6;  /* Wholesale Trade */
	if 52<=ind<=59 then ind=7;  /* Retail Trade */
	if 60<=ind<=67 then ind=8;  /* Finance, Insurance & Real Estate */
	if 70<=ind<=81 then ind=9;  /* Services */
	if 83<=ind<=89 then ind=9;  /* Services */
	if ind=82 then ind=10;  /* Education */
	if 90<=ind<=97 then ind=11;  /* Public Administration */
	if ind=99 | ind=0 | ind=. then ind=-9;  /* Unidentified */
run;
/***************************************/
/*  Compute # of firms each worker     */
/*  is employed at in an interval      */
/***************************************/
/*
proc sort data=fe_out; by id firmid; run;
data i_comp;
	set fe_out; by id firmid;
	if first.firmid then x=1;
	if first.id then ic_firm=0; ic_firm+x;
    if last.id;
	c_worker=ic_firm-1;
	keep id c_worker;
run;
data fe_out;
	merge fe_out i_comp;
	by id;
run;
*/
/***************************************/
/*  Compute # of comparison obs each   */
/*  firm has in an interval            */
/***************************************/
/*
proc sort data=fe_out; by firmid id; run;
data f_comp;
	set fe_out; by firmid id;
	if first.id then x=c_worker;
	if first.firmid then c_firm=0; c_firm+x;
    if last.firmid;
	keep firmid c_firm;
run;
data fe_out;
	merge fe_out f_comp;
	by firmid;
run;
proc datasets lib=work nolist;
	delete i_comp f_comp;
quit;
*/
/****************************************/
/*  Compute comparison firm groups      */
/****************************************/
/*
data fe_out;
	set fe_out;
	if 1<=c_firm<=5 then c_firm_g=1;
	if 6<=c_firm<=10 then c_firm_g=2;
	if 11<=c_firm<=20 then c_firm_g=3;
	if 21<=c_firm<=50 then c_firm_g=4;
	if 50<c_firm then c_firm_g=5;
run;
*/
/***************************************************/
/*  Compute worker-year weighted firm stats        */
/***************************************************/
proc means data=fe_out noprint nway;
	class firmid;
	var y pe fe xb r;
	output out=firm_stats1 mean(y)=fm_y mean(pe)=fm_pe mean(fe)=fm_fe mean(xb)=fm_xb mean(r)=fm_r N=n_wrk;
run;
/***************************************************/
/*  Compute firm weighted firm stats               */
/***************************************************/
proc sort data=fe_out; by firmid year id; run;
data firm_level;
	set fe_out;
	by firmid year;
	if first.year;
	keep firmid year firm_size;
run;
/* computes average firm size over the interval */
proc means data=firm_level noprint nway; 
	class firmid;
	var firm_size;
	output out=firm_stats2 mean(firm_size)=fm_size N=n_fper;
run;
/***************************************/
/*  Merge both firm stats together     */
/***************************************/
data firm_stats;
	merge firm_stats1 firm_stats2;
	by firmid;
	drop _type_ _freq_;
	if fm_size<=20 then fsize_g=1;    /* create firm-size group variable */
	if 20<fm_size<=100 then fsize_g=2;
	if 100<fm_size<=1000  then fsize_g=3;
	if 1000<fm_size<=10000  then fsize_g=4;
	if 10000<fm_size then fsize_g=5;
run;
proc datasets lib=work nolist;
	delete firm_stats1 firm_stats2;
quit;
data fe_out;  /* merge firm stats to main data */
	merge fe_out firm_stats;
	by firmid;
run;
/***********************************************/
/*  Macro that creates percentiles efficiently */
/*  	                                       */
/*  	- WORKER-YEAR WEIGHTED                 */
/***********************************************/
%macro rrank(ddata, vvar,percentile,groupvar,pername);

 proc sort data = &ddata out = &ddata._w; by &groupvar; run;
 
 proc sort data = &ddata._w;by &vvar; run;

 data &ddata._s; 
	set &ddata._w nobs=dsobst; 
	by &vvar;
    if _n_ = 1 then call symput ('dsobs',put(dsobst,13.));
    retain cntr 0 const1 cntrt;
    if first.&vvar then do; cntrt = 0; const1 = 0; end;
    cntr + 1;
    cntrt = cntrt + cntr;
    const1 + 1;
    if last.&vvar then do;
    mean_rank = cntrt / const1;
 output;
 end;
	keep &vvar mean_rank;
 run;

 proc sort data = &ddata ; by &vvar; run;

 data &ddata;
 	merge &ddata &ddata._s;
 	by &vvar;
 	rank=mean_rank/&dsobs.;
 	&vvar._&pername=ceil(rank*&percentile);
 	drop rank mean_rank;
	if &vvar._&pername=. then &vvar._&pername=0;
 run;
proc datasets lib=work nolist;
	delete &ddata._w &ddata._s;
quit;
%mend;
/***************************************/
/*  Create deciles for worker and      */
/*  firm fixed effects.                */
/*                                     */
/*  	- WORKER-YEAR WEIGHTED         */
/***************************************/
%rrank (fe_out,pe,10,id,dec);
%rrank (fe_out,fe,10,firmid,dec);
/****************************************/
/*  Create firm deviation variables    */
/****************************************/
/*
data fe_out;
	set fe_out;
	df_y=y-fm_y; df_pe=pe-fm_pe; df_fe=fe-fm_fe; df_xb=xb-fm_xb; df_r=r-fm_r; df_m=m-fm_m; df_e=e-fm_e;
run;
*/
/***********************************************/
/*  Create variables for all potential moments */
/***********************************************/
data fe_out;
	set fe_out;
	age_age=age*age;
	fsize_fsize=firm_size*firm_size; 
	fsize_y=firm_size*y; fsize_pe=firm_size*pe; fsize_fe=firm_size*fe; fsize_xb=firm_size*xb; fsize_r=firm_size*r; fsize_age=firm_size*age;
	y_y=y*y; pe_pe=pe*pe; fe_fe=fe*fe; xb_xb=xb*xb; r_r=r*r;
	y_pe=y*pe; y_fe=y*fe; y_xb=y*xb; y_r=y*r;
	pe_fe=pe*fe; pe_xb=pe*xb; pe_r=pe*r;
	fe_xb=fe*xb; fe_r=fe*r;
	xb_r=xb*r;
run;
/* drop some datasets */
proc datasets lib=work nolist;
	delete firmsizes firm_level firm_stats;
quit;
/*************************************************/
/*                                               */
/*  1.) STATS BY # OF JOB SWITCHES CATEGORY      */
/*                                               */
/*          - also by IND & FSIZE GRP            */
/*                                               */
/*************************************************/
/*
proc means data=fe_out noprint nway;
	class c_worker ind fsize_g;
	output out=swt_dist N=n_wkryears 
					mean(c_firm)=m_cfirm mean(n_jobs)=m_njobs
					mean(age)=m_age mean(age_age)=m_age_age
					mean(firm_size)=m_fsize mean(fsize_fsize)=m_fsize_fsize
					mean(y)=y mean(y)=m_y_y
					mean(pe)=m_pe mean(pe_pe)=m_pe_pe
					mean(fe)=m_fe mean(fe_fe)=m_fe_fe 
					mean(xb)=m_xb mean(xb_xb)=m_xb_xb 
					mean(r)=m_r mean(r_r)=m_r_r 
					mean(m)=m_m mean(m_m)=m_m_m 
					mean(e)=m_e mean(e_e)=m_e_e 
					mean(fm_pe)=m_fmpe mean(fm_pe)=m_fmpe_fmpe
					mean(pe_fe)=m_pe_fe;
run;
data swt_dist;
	retain interval;
	set swt_dist ;
	interval=&K;
	drop _type_ _freq_;
run;
*/
/*********************************************************/
/*                                                       */
/*  2.) STATS BY JOINT DIST OF WORKER & FIRM FE DECILES  */
/*                                                       */
/*          - also by IND & FSIZE GRP                    */
/*                                                       */
/*********************************************************/
proc means data=fe_out noprint nway;
	class pe_dec fe_dec ind fsize_g;
	output out=joint N=n_wkryears 
					mean(n_jobs)=m_njobs
					mean(age)=m_age mean(age_age)=m_age_age
					mean(firm_size)=m_fsize mean(fsize_fsize)=m_fsize_fsize
					mean(y)=y mean(y)=m_y_y
					mean(pe)=m_pe mean(pe_pe)=m_pe_pe
					mean(fe)=m_fe mean(fe_fe)=m_fe_fe 
					mean(xb)=m_xb mean(xb_xb)=m_xb_xb 
					mean(r)=m_r mean(r_r)=m_r_r 
					mean(pe_fe)=m_pe_fe;
run;
data joint;
	retain interval;
	set joint;
	interval=&K;
	drop _type_ _freq_;
run;
/*********************************************************/
/*                                                       */
/*  3.) WI-/BT-FIRM STATS BY IND, FSIZE GRP,             */
/*			& COMPARISON FIRM GRP                        */
/*                                                       */
/*********************************************************/
/*
proc means data=fe_out noprint nway;
	class ind fsize_g c_firm_g;
	output out=bt_wi_moments 
		N=n_wrkyear
		mean(c_worker)=m_cworker mean(c_firm)=m_cfirm mean(n_jobs)=m_njobs
		mean(y)=m_y mean(pe)=m_pe mean(fe)=m_fe mean(xb)=m_xb mean(r)=m_r mean(m)=m_m mean(e)=m_e
		mean(fm_y)=m_fmy mean(fm_pe)=m_fmpe mean(fm_fe)=m_fmfe mean(fm_xb)=m_fmxb mean(fm_r)=m_fmr mean(fm_m)=m_fmm mean(fm_e)=m_fme
		mean(df_y)=m_dfy mean(df_pe)=m_dfpe mean(df_fe)=m_dffe mean(df_xb)=m_dfxb mean(df_r)=m_dfr mean(df_m)=m_dfm mean(df_e)=m_dfe
		mean(y_y)=m_y_y mean(pe_pe)=m_pe_pe mean(fe_fe)=m_fe_fe mean(xb_xb)=m_xb_xb mean(r_r)=m_r_r mean(m_m)=m_m_m mean(e_e)=m_e_e
		mean(fmy_fmy)=m_fmy_fmy mean(fmpe_fmpe)=m_fmpe_fmpe mean(fmfe_fmfe)=m_fmfe_fmfe mean(fmxb_fmxb)=m_fmxb_fmxb mean(fmr_fmr)=m_fmr_fmr mean(fmm_fmm)=m_fmm_fmm mean(fme_fme)=m_fme_fme
		mean(dfy_dfy)=m_dfy_dfy mean(dfpe_dfpe)=m_dfpe_dfpe mean(dffe_dffe)=m_dffe_dffe mean(dfxb_dfxb)=m_dfxb_dfxb mean(dfr_dfr)=m_dfr_dfr mean(dfm_dfm)=m_dfm_dfm mean(dfe_dfe)=m_dfe_dfe
		mean(y_pe)=m_y_pe mean(y_fe)=m_y_fe mean(y_xb)=m_y_xb mean(y_r)=m_y_r mean(y_m)=m_y_m mean(y_e)=m_y_e
		mean(y_fmy)=m_y_fmy mean(y_fmpe)=m_y_fmpe mean(y_fmfe)=m_y_fmfe mean(y_fmxb)=m_y_fmxb mean(y_fmr)=m_y_fmr mean(y_fmm)=m_y_fmm mean(y_fme)=m_y_fme
		mean(y_dfy)=m_y_dfy mean(y_dfpe)=m_y_dfpe mean(y_dffe)=m_y_dffe mean(y_dfxb)=m_y_dfxb mean(y_dfr)=m_y_dfr mean(y_dfm)=m_y_dfm mean(y_dfe)=m_y_dfe
		mean(pe_fe)=m_pe_fe mean(pe_xb)=m_pe_xb mean(pe_r)=m_pe_r mean(pe_m)=m_pe_m mean(pe_e)=m_pe_e
		mean(pe_fmy)=m_pe_fmy mean(pe_fmpe)=m_pe_fmpe mean(pe_fmfe)=m_pe_fmfe mean(pe_fmxb)=m_pe_fmxb mean(pe_fmr)=m_pe_fmr mean(pe_fmm)=m_pe_fmm mean(pe_fme)=m_pe_fme
		mean(pe_dfy)=m_pe_dfy mean(pe_dfpe)=m_pe_dfpe mean(pe_dffe)=m_pe_dffe mean(pe_dfxb)=m_pe_dfxb mean(pe_dfr)=m_pe_dfr mean(pe_dfm)=m_pe_dfm mean(pe_dfe)=m_pe_dfe
		mean(fe_xb)=m_fe_xb mean(fe_r)=m_fe_r mean(fe_m)=m_fe_m mean(fe_e)=m_fe_e
		mean(fe_fmy)=m_fe_fmy mean(fe_fmpe)=m_fe_fmpe mean(fe_fmfe)=m_fe_fmfe mean(fe_fmxb)=m_fe_fmxb mean(fe_fmr)=m_fe_fmr mean(fe_fmm)=m_fe_fmm mean(fe_fme)=m_fe_fme
		mean(fe_dfy)=m_fe_dfy mean(fe_dfpe)=m_fe_dfpe mean(fe_dffe)=m_fe_dffe mean(fe_dfxb)=m_fe_dfxb mean(fe_dfr)=m_fe_dfr mean(fe_dfm)=m_fe_dfm mean(fe_dfe)=m_fe_dfe
		mean(xb_r)=m_xb_r mean(xb_m)=m_xb_m mean(xb_e)=m_xb_e
		mean(xb_fmy)=m_xb_fmy mean(xb_fmpe)=m_xb_fmpe mean(xb_fmfe)=m_xb_fmfe mean(xb_fmxb)=m_xb_fmxb mean(xb_fmr)=m_xb_fmr mean(xb_fmm)=m_xb_fmm mean(xb_fme)=m_xb_fme
		mean(xb_dfy)=m_xb_dfy mean(xb_dfpe)=m_xb_dfpe mean(xb_dffe)=m_xb_dffe mean(xb_dfxb)=m_xb_dfxb mean(xb_dfr)=m_xb_dfr mean(xb_dfm)=m_xb_dfm mean(xb_dfe)=m_xb_dfe
		mean(r_m)=m_r_m mean(r_e)=m_r_e
		mean(r_fmy)=m_r_fmy mean(r_fmpe)=m_r_fmpe mean(r_fmfe)=m_r_fmfe mean(r_fmxb)=m_r_fmxb mean(r_fmr)=m_r_fmr mean(r_fmm)=m_r_fmm mean(r_fme)=m_r_fme
		mean(r_dfy)=m_r_dfy mean(r_dfpe)=m_r_dfpe mean(r_dffe)=m_r_dffe mean(r_dfxb)=m_r_dfxb mean(r_dfr)=m_r_dfr mean(r_dfm)=m_r_dfm mean(r_dfe)=m_r_dfe
		mean(m_e)=m_m_e
		mean(m_fmy)=m_m_fmy mean(m_fmpe)=m_m_fmpe mean(m_fmfe)=m_m_fmfe mean(m_fmxb)=m_m_fmxb mean(m_fmr)=m_m_fmr mean(m_fmm)=m_m_fmm mean(m_fme)=m_m_fme
		mean(m_dfy)=m_m_dfy mean(m_dfpe)=m_m_dfpe mean(m_dffe)=m_m_dffe mean(m_dfxb)=m_m_dfxb mean(m_dfr)=m_m_dfr mean(m_dfm)=m_m_dfm mean(m_dfe)=m_m_dfe
		mean(e_fmy)=m_e_fmy mean(e_fmpe)=m_e_fmpe mean(e_fmfe)=m_e_fmfe mean(e_fmxb)=m_e_fmxb mean(e_fmr)=m_e_fmr mean(e_fmm)=m_e_fmm mean(e_fme)=m_e_fme
		mean(e_dfy)=m_e_dfy mean(e_dfpe)=m_e_dfpe mean(e_dffe)=m_e_dffe mean(e_dfxb)=m_e_dfxb mean(e_dfr)=m_e_dfr mean(e_dfm)=m_e_dfm mean(e_dfe)=m_e_dfe
		mean(fmy_fmpe)=m_fmy_fmpe mean(fmy_fmfe)=m_fmy_fmfe mean(fmy_fmxb)=m_fmy_fmxb mean(fmy_fmr)=m_fmy_fmr mean(fmy_fmm)=m_fmy_fmm mean(fmy_fme)=m_fmy_fme
		mean(fmy_dfy)=m_fmy_dfy mean(fmy_dfpe)=m_fmy_dfpe mean(fmy_dffe)=m_fmy_dffe mean(fmy_dfxb)=m_fmy_dfxb mean(fmy_dfr)=m_fmy_dfr mean(fmy_dfm)=m_fmy_dfm mean(fmy_dfe)=m_fmy_dfe
		mean(fmpe_fmfe)=m_fmpe_fmfe mean(fmpe_fmxb)=m_fmpe_fmxb mean(fmpe_fmr)=m_fmpe_fmr mean(fmpe_fmm)=m_fmpe_fmm mean(fmpe_fme)=m_fmpe_fme
		mean(fmpe_dfy)=m_fmpe_dfy mean(fmpe_dfpe)=m_fmpe_dfpe mean(fmpe_dffe)=m_fmpe_dffe mean(fmpe_dfxb)=m_fmpe_dfxb mean(fmpe_dfr)=m_fmpe_dfr mean(fmpe_dfm)=m_fmpe_dfm mean(fmpe_dfe)=m_fmpe_dfe
		mean(fmfe_fmxb)=m_fmfe_fmxb mean(fmfe_fmr)=m_fmfe_fmr mean(fmfe_fmm)=m_fmfe_fmm mean(fmfe_fme)=m_fmfe_fme
		mean(fmfe_dfy)=m_fmfe_dfy mean(fmfe_dfpe)=m_fmfe_dfpe mean(fmfe_dffe)=m_fmfe_dffe mean(fmfe_dfxb)=m_fmfe_dfxb mean(fmfe_dfr)=m_fmfe_dfr mean(fmfe_dfm)=m_fmfe_dfm mean(fmfe_dfe)=m_fmfe_dfe
		mean(fmxb_fmr)=m_fmxb_fmr mean(fmxb_fmm)=m_fmxb_fmm mean(fmxb_fme)=m_fmxb_fme
		mean(fmxb_dfy)=m_fmxb_dfy mean(fmxb_dfpe)=m_fmxb_dfpe mean(fmxb_dffe)=m_fmxb_dffe mean(fmxb_dfxb)=m_fmxb_dfxb mean(fmxb_dfr)=m_fmxb_dfr mean(fmxb_dfm)=m_fmxb_dfm mean(fmxb_dfe)=m_fmxb_dfe
		mean(fmr_fmm)=m_fmr_fmm mean(fmr_fme)=m_fmr_fme
		mean(fmr_dfy)=m_fmr_dfy mean(fmr_dfpe)=m_fmr_dfpe mean(fmr_dffe)=m_fmr_dffe mean(fmr_dfxb)=m_fmr_dfxb mean(fmr_dfr)=m_fmr_dfr mean(fmr_dfm)=m_fmr_dfm mean(fmr_dfe)=m_fmr_dfe
		mean(fmm_fme)=m_fmm_fme
		mean(fmm_dfy)=m_fmm_dfy mean(fmm_dfpe)=m_fmm_dfpe mean(fmm_dffe)=m_fmm_dffe mean(fmm_dfxb)=m_fmm_dfxb mean(fmm_dfr)=m_fmm_dfr mean(fmm_dfm)=m_fmm_dfm mean(fmm_dfe)=m_fmm_dfe
		mean(fme_dfy)=m_fme_dfy mean(fme_dfpe)=m_fme_dfpe mean(fme_dffe)=m_fme_dffe mean(fme_dfxb)=m_fme_dfxb mean(fme_dfr)=m_fme_dfr mean(fme_dfm)=m_fme_dfm mean(fme_dfe)=m_fme_dfe
		mean(dfy_dfpe)=m_dfy_dfpe mean(dfy_dffe)=m_dfy_dffe mean(dfy_dfxb)=m_dfy_dfxb mean(dfy_dfr)=m_dfy_dfr mean(dfy_dfm)=m_dfy_dfm mean(dfy_dfe)=m_dfy_dfe
		mean(dfpe_dffe)=m_dfpe_dffe mean(dfpe_dfxb)=m_dfpe_dfxb mean(dfpe_dfr)=m_dfpe_dfr mean(dfpe_dfm)=m_dfpe_dfm mean(dfpe_dfe)=m_dfpe_dfe
		mean(dffe_dfxb)=m_dffe_dfxb mean(dffe_dfr)=m_dffe_dfr mean(dffe_dfm)=m_dffe_dfm mean(dffe_dfe)=m_dffe_dfe
		mean(dfxb_dfr)=m_dfxb_dfr mean(dfxb_dfm)=m_dfxb_dfm mean(dfxb_dfe)=m_dfxb_dfe
		mean(dfr_dfm)=m_dfr_dfm mean(dfr_dfe)=m_dfr_dfe
		mean(dfm_dfe)=m_dfm_dfe
		mean(age)=m_age mean(age_age)=m_age_age
		mean(firm_size)=m_fsize mean(fsize_fsize)=m_fsize_fsize
		mean(fsize_y)=m_fsize_y mean(fsize_pe)=m_fsize_pe mean(fsize_fe)=m_fsize_fe mean(fsize_xb)=m_fsize_xb mean(fsize_r)=m_fsize_r mean(fsize_m)=m_fsize_m mean(fsize_e)=m_fsize_e;
run;
*/
/* format */
/*
data bt_wi_moments;
	retain interval;
	set bt_wi_moments;
	interval=&K;
	drop _freq_ _type_;
run;
*/
/* EXPORT */
/*
proc export data=swt_dist
	outfile="&outputdir.\&sample._SwtDist_&spec._&labsuprest._&minwage._&earn._&gender._&K..csv"
	DBMS=CSV
	replace;
	putnames=yes;
run;
*/
proc export data=joint
	outfile="&outputdir.\&sample._JointDist_&spec._&labsuprest._&minwage._&earn._&gender._&K..csv"
	DBMS=CSV
	replace;
	putnames=yes;
run;
/*
proc export data=bt_wi_moments
	outfile="&outputdir.\&sample._BT_WI_Firm_Mmnts_&spec._&labsuprest._&minwage._&earn._&gender._&K..csv"
	DBMS=CSV
	replace;
run;
*/
%end;
%mend;
%feanalysis(&earnval,&genderval,&minwageval);

/*********************************************************/
/*                                                       */
/*  COMBINE OUTPUT FROM ALL INTERVALS                    */
/*                                                       */
/*  (IF THE PROGRAM COMPLETES THIS STEP, THEN            */
/*   THE INTERVAL-SPECIFIC RESULTS CAN BE DELETED)       */
/*                                                       */
/*********************************************************/
libname firmdata "&firmdir.";
%macro combine(earn,gender,minwage);
%do K=1 %to 5;   
/*
proc import datafile="&outputdir.\&sample._SwtDist_&spec._&labsuprest._&minwage._&earn._&gender._&K..csv" 
	out=swt_dist_&K
	DBMS=CSV
	replace;
	getnames=yes;
run;
*/
proc import datafile="&outputdir.\&sample._JointDist_&spec._&labsuprest._&minwage._&earn._&gender._&K..csv"
	out=joint_&K
	DBMS=CSV
	replace;
	getnames=yes;
run;
/*
proc import datafile="&outputdir.\&sample._BT_WI_Firm_Mmnts_&spec._&labsuprest._&minwage._&earn._&gender._&K..csv"
	out=bt_wi_moments_&K
	DBMS=CSV
	replace;
	getnames=yes;
run;
*/
%end;
/*
data swt_dist;
	set %do K=1 %to 5; swt_dist_&K %end;;
run;
*/
data joint;
	set %do K=1 %to 5; joint_&K %end;;
run;
/*
data bt_wi_moments;
	set %do K=1 %to 5; bt_wi_moments_&K %end;;
run;
*/
/* EXPORT */
/*
proc export data=swt_dist
	outfile="&outputdir.\&sample._SwtDist_&spec._&labsuprest._&minwage._&earn._&gender..csv"
	DBMS=CSV
	replace;
	putnames=yes;
run;
*/
proc export data=joint
	outfile="&outputdir.\&sample._JointDist_&spec._&labsuprest._&minwage._&earn._&gender..csv"
	DBMS=CSV
	replace;
	putnames=yes;
run;
/*
proc export data=bt_wi_moments
	outfile="&outputdir.\&sample._BT_WI_Firm_Mmnts_&spec._&labsuprest._&minwage._&earn._&gender..csv"
	DBMS=CSV
	replace;
run;
*/
%mend;
%combine(&earnval,&genderval,&minwageval);
